
quit;
Libname sortie 'C:\My Documents\VAR\Econometrics\Data\';

data Temp1;
set sortie.Countriesq_clean;

Rzero=0.01; /* Annual riskless rate */
if year>1998.4 then delete;   /* selects estimation interval */

if country ^='Belgium' then delete; /* Selects country to be studied */
deltaNO=dif(NO);

run;


proc standard data=Temp1 mean=0 out=Temp2;

/* de-means the variables listed after the var statement -and these variables only */

var CA1 deltaNO;

run;

 data Temp3;
set Temp2;

CA1_1=lag(CA1);
CA1_2=lag2(CA1);
CA1_3=lag3(CA1);
CA1_4=lag4(CA1);
CA1_5=lag5(CA1);
CA1_6=lag6(CA1);
CA1_7=lag7(CA1);
CA1_8=lag5(CA1);
CA1_9=lag6(CA1);
CA1_10=lag7(CA1);


deltaNO_1=lag(deltaNO);
deltaNO_2=lag2(deltaNO);
deltaNO_3=lag3(deltaNO);
deltaNO_4=lag4(deltaNO);
deltaNO_5=lag5(deltaNO);
deltaNO_6=lag6(deltaNO);
deltaNO_7=lag7(deltaNO);
deltaNO_8=lag5(deltaNO);
deltaNO_9=lag6(deltaNO);
deltaNO_10=lag7(deltaNO);



run;



proc varmax data=Temp2 outest=est outcov; 

model  CA1 deltaNO/method=ls noint noprint p=3 ;

run;



proc iml;
/* ----------------PARAMETERS------------ */
iterations=10000;
Rzerobis=0.01;
/* ----------------PARAMETERS------------ */
use est;
read all var _num_ into y where(type="EST"); 
nombrelag=3; 
nombrevar=2; 
nombrecoef=nombrelag*nombrevar*nombrevar;
y=y[1:nombrevar,1:nombrevar*nombrelag];
Aint1=I(nombrevar*(nombrelag-1));
Aint2=j(nombrevar*(nombrelag-1),nombrevar,0);
Aint=Aint1||Aint2; 
print y Aint;

A=y`||Aint`;

A=A`;



disc=1/(1+RZerobis);

h=J(1,ncol(y),0);
h[1,2]=-1; 
h=h*disc;

Id=I(nombrevar*nombrelag);
inverse=inv(Id-disc*A);

K=h*A*inverse;
print K;


use Temp3;
/* ----------------------PARAMETERS------------------ */
read all var  {CA1 deltaNO CA1_1 deltaNO_1 CA1_2 deltaNO_2} into g where (CA1^=. & deltaNO^=. & CA1_1^=. & deltaNO_1^=. & CA1_2^=. & deltaNO_2^=.);  
/* ----------------------PARAMETERS------------------ */
/*    {CA1 deltaNO CA1_1 deltaNO_1 CA1_2 deltaNO_2 CA1_3 deltaNO_3 CA1_4 deltaNO_4 CA1_5 deltaNO_5  CA1_6 deltaNO_6 CA1_7 deltaNO_7 CA1_8 deltaNO_8} into g where (CA1^=. & deltaNO^=. & CA1_1^=. & deltaNO_1^=. & CA1_2^=. & deltaNO_2^=.& CA1_3^=. & deltaNO_3^=.& CA1_4^=. & deltaNO_4^=. & CA1_5^=. & deltaNO_5^=. & CA1_6^=. & deltaNO_6^=. & CA1_7^=. & deltaNO_7^=. & CA1_8^=. & deltaNO_8^=.);    */

CAhat=K*G`;

h=J(1,ncol(G),0);
h[1,1]=1;

CA1=h*G`; 

un=J(ncol(Cahat),1,1); 
Cahatmean=(Cahat*un)/ncol(CAhat); /
Cahatvar=(CAhat*CAhat` -ncol(CAhat)*Cahatmean*Cahatmean)/(ncol(CAhat)-1);

CA1mean=(CA1*un)/ncol(CA1); /* moyenne de Cahat*/
CA1var=(CA1*CA1` -ncol(CA1)*CA1mean*CA1mean)/(ncol(CA1)-1);

CAcovar=(CAhat*CA1` -ncol(CAhat)*Cahatmean*CA1mean)/(ncol(CAhat)-1); 
CAcorr=CAcovar/sqrt(Cahatvar*CA1var);

print CA1mean Cahatmean Cahatvar CA1var CAcorr;

use Temp3;
/* ----------------------PARAMETERS------------------ */
read all var {year} into year where (CA1^=. & deltaNO^=.  & CA1_1^=. & deltaNO_1^=. & CA1_2^=. & deltaNO_2^=.); /* recupere la variable year */
/* ----------------------PARAMETERS------------------ */

create dessine var{CA1 CAhat year};
append;
close dessine;
sort dessine out=dessine by year;



use est;
read all var _num_ into varcovar where(type="COV"); 
varcovar=varcovar[1:nombrecoef,1+nombrevar*nombrelag:nombrecoef+nombrevar*nombrelag];


print varcovar; 

grandvect=(y[1:nombrevar,1])`; 

do i=2 to nombrevar*nombrelag;
		intermed =(y[1:nombrevar,i])`;
		grandvect=grandvect||intermed;
   end;

print grandvect;


grandK=j(iterations,ncol(K),0); /* initialisation*/
print(ncol(CAhat));
grandCAsim=j(iterations,ncol(CAhat),0);  /* initialisation*/


print varcovar;
h=J(1,ncol(y),0);
h[1,2]=-1; 
h=h*disc;

do iter=1 to iterations; 


grandvectsimul=(grandvect)`+ root(varcovar)*normal(j(nombrecoef,1,iter)); /* generates a multivariate normal vector of mean and variance as indicated */

ysimul=j(nombrevar,nombrevar*nombrelag,0); /

   do i=1 to nombrevar*nombrelag;
		ysimul[1:nombrevar,i]=grandvectsimul[1+(i-1)*nombrevar:i*nombrevar];
	end;



/* construction of the A matrix on stacked form */
A=ysimul`||Aint`;
A=A`;

inverse=inv(Id-disc*A);
Ksimul=h*A*inverse;
CAhatsim=Ksimul*G`;

Cahatsimmean=(Cahatsim*un)/ncol(CAhatsim); 
grandCAsim[iter,1:ncol(CAhat)]=Cahatsim;

end;


grandCAsimmean=(grandCAsim)`*j(iterations,1,1)/iterations;

print grandCAsimmean; 




create pourcentage from grandCAsim;
append from grandCAsim;
close pourcentage;

quit;




proc univariate data=pourcentage;  
/* ----------------------PARAMETERS------------------ */
var col1-col73; 
output out=location
 			pctlpts=2.5, 97.5
			pctlpre=col1-col73 
/* ----------------------PARAMETERS------------------ */
			pctlname=bas haut;

run;


proc iml; 
use location;
/* ----------------------PARAMETERS------------------ */
nombrevar=2; /* idem*/
/* ----------------------PARAMETERS------------------ */

read all var _num_ into y; 
print y;
CAbas=J(ncol(y)/2,1,0);
CAhaut=J(ncol(y)/2,1,0);

do i=1 to ncol(y)/2;
	CAbas[i,1]=y[1,2*i-1];
	CAhaut[i,1]=y[1,2*i];
end;

print CAbas CAhaut;

create hautbas var{CAbas CAhaut};
append;
close hautbas;


quit;




data dessine4;
   merge dessine hautbas;
  run;




symbol1 color=blue 
        interpol=join
        value=dot 
        height=0.5; 

symbol2 color=red 
        interpol=join
        value=triangle
        height=0.5; 


symbol3 color=green 
        interpol=join
        value=dot
		line=3
        height=0.1; 


symbol4 color=green 
		 interpol=join
        value=dot
		line=3
        height=0.1; 


proc gplot data=dessine4;
plot CA1*year Cahat*year CAbas*year CAhaut*year/ overlay haxis=1980 to 2002 by 5 
					vref=0; 
run;

proc gplot data=dessine4;
plot CA1*year  Cahat*year/ overlay haxis=1980 to 2002 by 5 

					vref=0; 
run;
proc gplot data=dessine4;
plot CA1*year CAbas*year CAhaut*year / overlay haxis=1980 to 2002 by 5 
					vref=0; 
run;
proc gplot data=dessine4;
plot Cahat*year CAbas*year CAhaut*year/ overlay haxis=1980 to 2002 by 5 
					vref=0; 
run;
quit;
      


